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Recent discovery of kilohertz quasi-periodic brightness oscillations of low mass X-ray 
tj^l binaries (LMXBs) has attracted attention to highly relativistic periodic motion near 

^^ ■ accreting neutron stars. Most models proposed so far involve (almost) inertial motion 

in the vicinity of the stars' innermost stable circular orbits. In the present paper we 
study general-relativistic circular and eccentric orbits around spinning neutron stars 
assuming the orbits are slightly tilted with respect to the stars' equatorial planes. We 
develop analytical and numerical techniques for integrating bound timelike geodesies 
^ . in fully relativistic neutron star spacetimes obtained by modern numerical codes. We 

f~^ ' use equations of state of neutron star matter that span a broad range of stiffness, 

V/^ \ while the explored range of masses (M > I.TMq) and spin frequencies [v^ < 600 Hz) 

■t;:J" . is motivated by the observations of LMXBs. We investigate the general-relativistic 

0^ ' effects of periastron advance and nodal precession in the strong gravitational fields of 

^D , rotating neutron stars and compare quantitatively the associated orbital frequencies 

w/ ■ with the more readily obtainable frequencies of orbits around Kerr black holes on 

the one hand, and low-order post-Newtonian (PN) expansions, on the other. While 
C^ , Kerr results approximate the periastron advance frequency much better than the PN 

^ i' expressions, the retrograde torque caused by the high quadrupole moments of the 

rotating stars clearly favours the PN approximation in the case of nodal precession. 
The methods developed in the present paper are used in the companion paper to 
test recent hypotheses regarding the origin of quasi-periodic oscillations in accreting 



o 

j^ ', neutron star sources. 



^ 
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1 INTRODUCTION AND OVERVIEW the definition and properties of the Z and atoll sources see 

, , T^ -.^rriT^ Hasinger & van der Klis 1989) have been observed at lower 

Over the past tour years, observations by the Rossi XTE „ ■ i^tt ^ ^ r.r,TT rm r ■ ^ 

... . , ,r , ■ ■ /T /i-iT-r. N 1 trequencies 10 tlz;$ j^l ;$dOHz. inese trequencies also vary 

satellite oi low-mass X-ray binaries (LMXBs) nave pro- 



in apparently tight correlation with the kHz QPO frequen- 
cies. 



duced an abundance of data on the rapid variations of these 
sources' X-ray brightness (see van der Klis 1998, 2000, for 
an overview, and references therein). Most prominent in The fact that the kHz QPO frequencies lie in the range 
the sources' high-resolution power spectra, the kilohertz- of orbital frequencies of bound geodesies near neutron-star- 
range quasi-periodic oscillations (QPOs) have been detected mass objects [the Kepler frequency for the innermost stable 
in ~ 20 binaries with compact accretors identified as neu- circular orbit (ISCO) around a static, spherical star of mass 
tron stars. The kHz QPOs often occur in pairs at frequen- M is i/isco = 1100 Hz (2Mq/M)] motivated several theoreti- 
cies v\ and V2 > vi (with V2 ranging between 500 Hz and cal scenarios (see, e.g., Strohmayer et al. 1996, Miller, Lamb 
1200 Hz) that can vary in a single source by up to a fac- & Psaltis 1998, Stella & Vietri 1999) for the X-ray luminos- 
tor of ~ 2 (over as little as few minutes) in step with an ity modulation that all involve inertial or near-inertial mo- 
inferred change in the accretion rate. In addition to the tion of localised condensations — variously named 'clumps' 
kHz QPOs, the horizontal-branch oscillations (HBOs) of or 'blobs' — in the accretion fiows around neutron stars. For 
the high-luminosity (near Eddington) neutron-star Z sources instance, in a purely inertial (geodesic) orbital model, Stella 
and 'bumps' in the power spectra of the lower-luminosity & Vietri (1999; see also Karas 1999) have proposed for a 
(between 10~^Z/Edd and a few 10~^LEdd) atoll sources (for given source an identification of the variable QPO frequen- 
cies 1^2 and 1^1 with, respectively, the Kepler frequencies uk 
and the periastron advance frequencies i/p^ of a family of ec- 
* Also Department of Physics. centric orbits. The preferred plane of accretion (and, thus, of 
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the clumps' motion) is normally assumed to lie close to the 
rotation equator of the central neutron star, a natural con- 
sequence of accretion-driven neutron star spin-up (similar 
reasoning clearly favours prograde motion in the accretion 
flow). However, we need not exclude small deviations of the 
orbits from the equatorial plane, and so it has been specu- 
lated that a signal — e.g., the QPOs observed at frequency 
i/L — might be produced by the nodal precession of the or- 
bital planes at frequency u^p or its first overtone (see, e.g., 
SteUa & Vietri 1998 and Morsink & Stella 1999 for the nodal 
precession of circular orbits) . 

If any of the orbital models is confirmed, observation 
of orbital frequencies would provide a direct probe of the 
highly curved spacetimes around rapidly spinning neutron 
stars and an unprecedented test of the general theory of rel- 
ativity in the strong-field regime. The highly accurate mea- 
surements of the kHz QPO frequencies (e.g., the relative er- 
ror < 1% for Sco X-1; van der Klis et al. 1997) impose, how- 
ever, stringent requirements on the QPO models; they also 
make it necessary to know the theoretical orbital frequencies 
with an accuracy that may render post-Newtonian approx- 
imations inadequate. Similarly, the spacetimes around neu- 
tron stars can differ from the Kerr black-hole spacetimes to 
such a degree as to invalidate the use of Kerr-based orbital 
frequencies. 

It is our purpose in the present and a companion paper 
to investigate in some detail orbital motion near the equa- 
torial planes of neutron stars and to test carefully the agree- 
ment between the frequencies of bound orbits with the QPO 
frequencies observed so far in neutron-star LMXBs. Specific 
QPO models will be compared with the existing data in the 
companion paper (Markovic & Lamb 2000). In the present 
paper we study the motion of test particles on orbits of ar- 
bitrary eccentricity around rapidly spinning neutron stars. 
We pay particular attention to the general-relativistic effects 
of periastron advance and nodal precession in the neutron 
stars' strong gravitational fields, and compare quantitatively 
the orbital frequencies associated with these effects with the 
more readily obtainable frequencies of orbits around Kerr 
black holes on the one hand, and low-order post-Newtonian 
(PN) expansions, on the other. As the basis of our study, 
we set out to develop analytical and numerical techniques 
for a sufficiently accurate computation of the frequencies 
of orbits in fully general relativistic numerical neutron star 
spacetimes. Modern numerical codes (Komatsu, Eriguchi & 
Hachisu 1989a,b; Cook, Shapiro & Teukolsky 1992, 1994a,b; 
Stergioulas & Friedman 1995; Nozawa et al. 1998) allow 
an accurate computation of equilibrium uniformly spinning 
neutron star models starting from an arbitrary tabulated 
equation of state (EOS), and it is one variant (Stergioulas 
& Friedman 1995) that we use to obtain external neutron 
star spacetimes in which we then integrate the equations of 
motion for bound geodesies. 

The numerical neutron star models discussed in this pa- 
per are all based on three equations of state that sample a 
broad range of stiffness, from the moderately stiff EOS C 
(Bethe & Johnson 1974; see also the early revi ew by Ar- 
nett fc Bowers, 1977 ) to the very stiff EOS L ( |Pandhari 



by the strong short-range repulsion due to the w-meson ex- 
change. EOS C allows the maximum mass of 1.86Af0 for 
non-rotating stars. On the other hand, EOS L, based on 
coupling of (non-relativistic) nuclear matter to a mean scalar 
field, provides for a much higher maximum mass, 2.71Mq, 
for non-rotating neutron stars and, accordingly, larger equa- 
torial radii. For the third EOS of intermediate stiffness we 



choose EOS UU (Wiringa, Fiks fc Fabrocini 198^), which. 



besides including the nucleon-nucleon potential (specifically, 
the Urbana U14), achieves the relatively greater stiffness at 
high density due to the Urbana VII three-nucleon poten- 
tial (Schiavilla, Pandharipande, fc Wiringa 1986) . The max- 
imum non-rotating NS mass for this EOS is 2.2OM0. [The 
more recent EOS A18-I-UIX -f- Svt of Akmal, Pandhari- 
pande and Ravenhall (1998; see also Pandharipande, Akmal 
& Ravenhall, 1998) is based on the modern Argonne vis 
two-nucleon potential and the Urbana IX three-nucleon po- 
tential and allows for the dependence of the interaction be- 
tween two nucleons on their total centre-of-mass momentum; 
its stiffness is similar to that of EOS UU and the maximum 
non-rotating stellar mass is 2.23AIq.] These three equations 
of state span the entire range of EOS stiffness that can pro- 
vide the neutron star masses (1.7 < M < 2.2; from now 
on all masses will be given in units of solar mass) and radii 
require d by the QPO models d iscussed in the companion 
paper (Markovic & Lamb 200C). 



pande fc Smith 1975). The classic EOS C describes cold nu- 
cleon (n and p) + hyperon (A, E ''^, A ''^ and A^""") matter 
with the nucleon-nucleon potential represented by a set of 
Yukawa functions. The high-density behaviour is dominated 



In Section 2, we derive the equations of motion for a 
general bound geodesic near the equatorial plane of the sta- 
tionary spacetime around a neutron star. The equations of 
motion involve the values of four radius-dependent poten- 
tials (provided as a tabulated output from a NS-model nu- 
merical code) in the equatorial plane, as well as the poten- 
tials' second derivatives with respect to the vertical coordi- 
nate. The presence of two integrals of motion, the energy E 
and the vertical component, L, of the angular momentum, 
allows integration of the radial and azimuthal equations up 
to quadratures (we specifically solve for prograde orbits). 
For orbits (slightly) tilted relative to the equatorial plane, 
the lack of a third integral of motion leaves one with two 
first-order Hamilton's equations for the vertical 0— motion. 
The precession of the orbital plane is then most conveniently 
analysed by means of Floquet's method. In the balance of 
Section 2 we discuss similar motion in the Kerr black hole 
spacetime, where the presence of a third integral of motion, 
the Carter constant, simplifies the treatment considerably. 
Finally, we briefly discuss low-order post-Newtonian expres- 
sions for the periastron advance and nodal precession fre- 
quencies. 

In Section 3, we illustrate, by means of several exam- 
ples, the salient features of the orbits around neutron stars, 
devoting most of our discussion to the dependence of the 
orbital frequencies, vk, I'r = vk — J^pa and i'npi on the radii 
of circular orbits or the periastra and apastra of eccentric 
orbits (unless specifled otherwise, all radii will be given in 
units GM/c?). We pay particular attention to the manner 
in which the frequency profiles are affected by equations of 
state of various stiffness. For instance, the stiffer EOS yield 
neutron stars of greater extent, which, in turn, can exclude 
highly eccentric orbits with small periastra and low Vr. The 
larger stellar extent (and, accordingly, greater magnitude 
of the quadrupole moment qG^ M'^ / c^) also reduces signif- 
icantly the nodal precession frequencies for stars of given 
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mass M and angular momentum jGM^ /c (this effect was 
previously explored for circular orbits by Morsink & Stella 
1999). We also assess the accuracy of approximating the 
orbital frequencies in neutron-star spacetimes by the more 
readily computed frequencies of orbits in Kerr spacetimes 
of the same mass M and angular momentum jGM^ /c, or 
by the frequencies obtained from low-order post-Newtonian 
expressions (including the relatively large mass quadrupole 
terms characteristic of neutron stars). The results of this 
paper are summarised and discussed in Section 4. 

Equipped with accurately computed frequencies of or- 
bits in neutron star spacetimes, we proceed in the companion 
paper (Markovic & Lamb 200C) to investigate the possibil- 



ity of their matching the QPO frequencies of neutron-star 
LMXBs. 



BOUND GEODESICS AROUND SPINNING 
COMPACT STARS 



2.1 Neutron stars 



The code of Stergioulas (Stcrgioulas & Friedman 1991: 



yields the metric around a neutron star in the form isotropic 
in the f — 9 sector 



ds 



„P+7 



dt^ + e^'Pf^ sin^ e{d(j> - ujdtf 



+e''°'(df^ + f^de^) 



(1) 



where ^{f,6), p{f,6), a{f,6), Lu{f,9) are the four poten- 
tials computed numerically. The integration of orbital equa- 
tions, described below, requires the four potentials to be 
given as smooth functions of f . We achieve sufficient smooth- 
ness through a combination of spline and high-order Cheby- 
shev polynomial representations (see, e.g.. Press et al. 1992). 
Throughout the paper, r will denote only the coordinate- 
independent circumferential radius. For metric (|l[), r = 
exp[(7 - p)/2)]f. 

The metric components' independence of t and <j) im- 
plies that the components ut = —E and u^ = L of the 
4-velocity of a test particle are conserved. It immediately 
follows that 

— - --(''+^'(£-a;L) 

^ "-^ ■ -'•''+''\E-luL). (2) 



dr 
d4 
d7 



r^ sin y 

Assume for the moment a purely equatorial motion, 
6 = 7r/2. The normalisation u'^m^ = — 1 leads then to the 
equation for radial motion 

- F(f) 



m 



-(P+7) 



{E- 



- ^Lf - L^e"-'' 



(3) 



For given apastron radius fa and periastron radius fp, 
a bound orbit exists if one can find L {L > for prograde 
orbits) such that F^, = Fp = (subscripts denote the values 
of f at which the quantities are evaluated), i.e., such that 
solve the equations 



E = UJ^L+ (AaL^ 



Ba 



,1/2 



JpL + {ApL^ + Bp 



,1/2 



(4) 



where A{r) = exp[2p(r)]/r^ and B{r) = exp[p(r) + 7(j')]. 
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Figure 2. The drift of ascending nodal points for eccentric or- 
bits of common periastron rp = 4.7 (the same NS as in Fig. 1; 
throughout the paper r denotes the circumferential radius) and 
apastra ra = 8.0 (uk = 1400 Hz, u^-p = 38 Hz), 9.5 (1120, 26), 
11.3 (964, 21), 13.4 (821, 17), 16.0 (692, 14), 19.2 (577, 11) and 
23.0 (475, 8.6). 

The required solution exists only if fp > fmb, where fmb 
is the periastron of the marginally bound {E = 1) orbit. 
Generally, for L and E corresponding to a bound orbit, there 
also is a third root of the equation f (fa) ~ 0, closer to the 
origin, fa < fp. For the marginally bound orbit fp = fa, 
f a = C)0 [-F'(fp) = 0, F"(fp) > 0; the prime denotes the 
derivative with respect to f]. 

For a given f p < f isco , where f isco is the radius of the in- 
nermost stable circular orbit, there is a minimum apastron, 
fa,miii(fp) > fisco below which there exist no bound (pro- 
grade) orbits. For fp > fisco we naturally have fa,min = fp. 

Once L and E are found for given fp and fa, it is 
straightforward to obtain the period (in the distant ob- 
server's time i) of radial motion 



— = 2 



E-ujL 



-(p+7)^^ 



(5) 



as well as the frequency of azimuthal motion, i.e., the Kepler 
frequency 



n 
where 

A<^ = 



(6) 



— dr 
r 



'-e"-'' + ioe-'-''+'<'>{E-LuL) 
The periastron advance frequency is then simply 



dr 



J'PA ^ V\i — Vr 



(7) 



(8) 



In the more general case of orbits out of the equatorial 
plane, the motion is governed by the Hamiltonian (see, e.g.. 
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Figure 1. Plots of cost/) (dotted lines) and •ip/4'o (solid lines) for three eccentric orbits of common periastron rp = 4.7 around a 
neutron star of mass M = 1.8, EOS UU, u^ = 500 Hz (j = 0.20). The apastra and the orbital frequencies are Ta = 8.0 (i'k = 1400 Hz, 
i/„p = 38 Hz), 13.4 (821, 17) and 23.0 (475, 8.6). 



Misner, Thorne & Wheeler 1973) 



H = -ie-<-''+-'^{E-u;Lf + le''-'<_^ 



, 1 -2a f 2 , 1 



2 



(9) 



For small displacements, ij) = n/2 — 9 (p^ = —pe), from 
the equatorial plane, we can expand the four potentials, 
p{f,ilj) « p(f,0) + p(f)V'V2 = Poir) + p{f)ip^/2, etc., and 
thus arrive at the 'harmonic oscillator' form quadratic in 
both %{) and p^ 



H = Ho{f,Pf) + \^\f)^^ + 



2/(f 



2 
■Pi,, 



I{f) = e^'^^'f^ 






and 




^"(f) = e-'^+^o) 


{E^ 


-oL)^(f + i-^) 




-^{E - uqL)Lu} 



(10) 



(11) 



+e"^(l + f-I+.)+, 



(12) 



Since the radial motion is (to the lowest order in ^l>) given 
by equation (H) [and thus Ha{f,pr) = —1/2], we have used 
this fact to substitute pr — exp(2Qo)df/dr in the expression 
for n'^{f). 

The motion in the i/; direction can now be computed 



from Hamilton's equations, given here in the matrix form. 



d_ 
d7 



i, 
Pi, 



Pi, 



Ai2 = l//(r), 



A2 



-n'if) 



All 



0, (13) 



where f(r) is a solution of equation (Bl). 

The secular advance of the nodal line (i.e., the nodal 
precession) is a straightforward consequence of the failure 
of a particle to return to the initial ^ after a full revolu- 
tion around the star. For example, in Fig. 1 we plot cos (j> 
and il^/^o for three orbits of common periastron (circum- 
ferential) radius rp = 4.7, which is smaller than risco = 5.4 
for the neutron star of mass M = 1.8, EOS UU and spin 
frequency v^ = 500 Hz (j — 0.20). Notice the gradually in- 
creasing phase lag of the meridional motion relative to the 
azimuthal motion. This leads to a secular advance of the 
ascending node (see Fig. 2). 

For eccentric orbits, the rate of the advance varies some- 
what over a few orbital periods i^^^'- since the radial fre- 
quency Ur is much smaller than vk for orbits close to risco 
(e.g., in the above example for ra = 8.0, Vr = 223 Hz = 
i'k/6.3), a particle performs several revolutions (see the plot 
of cos (f) in Fig. 1) while lingering close to the periastron, and 
then moves closer to the apastron where both the rate of rev- 
olution and — in particular — the effects of the star's spin 
and quadrupole deformation (oblateness) on the orbital mo- 
tion are much weaker. However, when averaged over several 
orbital periods, the nodal precession rate is constant, and 
this averaged value is what we denote by v^p . 

Equations (|l3|) can be solved by direct integration over 



© 0000 RAS, MNRAS 000, 000-000 



Bound near- equatorial orbits around neutron stars 5 



a large number of (proper time) periods Tr of radial motion 
to obtain the nodal drift with sufficient accuracy, as was 
done for Fig. 2. More economically, the same result can be 
obtained by Floquet's method (see, e.g., Arnol'd 1992 or 
Hochstadt 1975), as follows: We define the matrix 



* 



Pi P2 



(14) 



whose components are the values of the latitude ^p and 
the associated momenta p^ at the end of a period of ra- 
dial motion, tpi = ip (tt) , pi = p^ (r^) for initial conditions 
V'(O) = 1, p^(0) = 0, and V'2 = ip (rr), P2 = Pi, (rr) for ini- 
tial conditions ^(0) = 0, p^(0) = 1. Since the product of the 
eigenvalues of the matrix ^, /ii/i2 = exp (/J'^ |Ajdrj — 1, 
we obtain 



^1 = M2 



g'"^/".^ CQs{v^/v^) = \ ^ , 



(15) 



where one must include (by an appropriate addition of n27r) 
the total phase change in tp accumulated over Tr when calcu- 
lating v^/vr from equation (|l5[). The nodal precession fre- 
quency i^p^ is then simply 



uk~ Vi;. 



(16) 



The constants of motion L and E for almost equatorial 
circular orbits can be found at all f > fisco from the equa- 
tions F{f) = F'{f) — 0. It is then straightforward to obtain 
the Kepler frequency [see equations (Q)] 

and the radial epicyclic frequency 



2.Ur = tOm^e^-^ 



(18) 



E-ujL 

Using equations (h_l|) and (O) we find the meridional fre- 
quency 



27r//^ 



■y + p — a 



(19) 



{E -uL)f 
and thus the nodal precession frequency u^^ = vk — v-iii- 

2.2 Kerr black holes 

Expressed in the Boyer-Lindquist coordinates, the motion 
of a massive particle near the equatorial plane of a Kerr 
black hole is governed by the following equations (see, e.g., 
Chandrasekhar 1983) 

.4X2 

r r — TZ, 

-4 ;2 r\ 

r ip —IS, 



L[2jEr + L{f''-2f)] 
l[i5E^-2,Lf], 



(20) 



where 

A = f^-2f + j\ Y?={r^+ff -fA (21) 

and 

7^(f ) = - (l - E^) f* + 2f'^ - \j^ (l - E^) + L^] f^ 

+2{L-jEfr, 
© 0000 RAS, MNRAS 000, 000-000 



eW = [L^+f{l-E')]{^^o-^')- (22) 

In the above equations we have set t/j = everywhere except 
for the terms of order ip'^ in the equation for the ?/j-motion. 
The amplitude ipo of the angular displacement i/) from the 
equatorial plane is related to the Carter constant. 

For bound orbits {E < 1), the equation TZ{r) = has 
three positive roots fa > rp > fa, and the first of equations 
(bol) takes on the form 

f'^f^ = - (l - E^) {r - fa)(f - fp)(f - fs). (23) 

Taking fa and fp as independent variables and defining 
yl = f a -I- f p and B = fafp, we obtain from the coefficients 
of equation (GSf) the following relations 

2 



•2 , 
J + 



l-£'2 
l-£2 

(L-jE)^ 

l-£2 



= A + fs, 
= B + Ars, 
= Bfi- 



(24) 



Eliminating E and L from equations (|24) we obtain a 
second-order equation for r^ 

Difi + Difa + Do = 0, (25) 

where 

Do = {2B~fAY 

Di = 2(2B-fA){2A-B-f)-4:f{A-2)B 

Di = {2A-B-fY~4fB. (26) 

Once fa is found for given fp and fa, the angular momentum 
and energy are obtained straightforwardly 

,2 2 



f a + f p + f 3 

2 



(fafp + fpfs + fafa - f) , 



E^ = I -. (27) 

ra + T-p + r-3 

Similarly to the neutron stars discussed above, for fp > 
fmb, where 



T-mb 



=2-j+2^T^ 



J 



(28) 



is the Boyer-Lindquist periastron of the marginally bound 
{E = 1) prograde orbit, eccentric orbits exist only if fa > 
fa.min where fa,min(fp) > fisco and 

fisco = 3+Z2-[(3-Zl)(3 + Zl+2Z2)]'^^ 



zi^i+{i-ff'[{i+jr'+{i-j)'^'] 



{3f + Z 



2\l/2 



(29) 



is the radius of the innermost stable circular orbit. Again, 

^a,min — ^p U /"p ^ ?^isco. 

The two positive roots (if they exist, i.e., if there are 
bound orbits with the given apastron fa and periastron fp) 
of equation (M) correspond to prograde and retrograde or- 
bits. Restricting ourselves for the moment to weakly bound 
orbits at large radii 



fa,fp > 1; 



1 + E: 



u = i-E<s:i, 



(30) 



around almost spherical, j <^ 1, black holes, we obtain from 
equations 
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"^'-^ B-2A + 2jL ■ ^^^' 

We see that f3{L > 0) < r3(L < 0) and that the pro- 
grade orbits are more tightly bound than the retrograde 
ones, U{L > 0) > U{L < 0). The prograde orbit thus corre- 
sponds to the smaUer of the roots of equation (Ea). 
The period of the radial motion is 



— = 2 



—dr 



ET? - 2jLf 



A^f(l-S2)(r-f3) 



dx, 



(32) 



where we have introduced a convenient integration parame- 
ter X such that 



r = 1{1 + ecosx), I 



ra + rp 



(33) 



2 ra -I- Tp 

Using the general relation (n) to find the Kepler frequency, 
we integrate the third of equations (EOh and thus obtain 



A(?i = 



— df 

r 



2jEr + i (f ^ - 2f ) 



--dx- 



(34) 



Finally, the second of equations (M) and the second of equa- 
tions d23) give the motion in the i/i-direction 



(35) 



tP = i>o cos (y/ L'^ + f {1 - E^) / Ij^dA . 

The failure of a particle to return to the initial i/) after a full 
period of azimuthal motion will lead to a secular precession 
of the orbital plane (i.e., the nodal precession) at frequency 
([16|) where 



v.^, 



Ur-VL-'+fil-E^) [ ' l^ 
^f^{l-E^){f-h) 



dr 



(36) 



For a circular orbit at Boyer-Lindquist radius f near the 
equatorial plane, the Kepler frequency for prograde orbits is 
(see, e.g., Chandrasekhar 1983) 



2'm'K 



p/2 (1+ jf-3/2-) 



(37) 



(to obtain expressions for retrograde orbits one should set 
j -^ —j)- At all f > fisco a slight perturbation of the orbital 
motion in the radial direction will lead to radial oscillations 
at the radial epicyclic frequency 

-2\l/2 



Ur ~ \i^k\ (i — Gf + 8jf 



-3/2 



o .2 .- 

3j r 



(38) 



On the other hand, a particle in a circular orbit tilted 
slightly with respect to the equatorial plane will undergo 
meridional oscillations at frequency 



leading to the nodal precession at frequency (uMj. 



(39) 



At large radii (r^lO), the asymptotic expressions 

2jf~^ lead to the simple 



2niy-p^ ^ 3f 
scaling 



-5/2 



and 27r!/,, 



"(?.\"'.."' 



""•(5) 



(40) 



At smaller radii, the slope d In i/^^p /d In i/p^ steadily de- 
creases with increasing frequencies and drops below unity 
near the innermost stable circular orbit. 



2.3 Spinning stars in post-Newtonian 
approximat ion 

At a sufficiently large distance from a spinning star of angu- 
lar momentum jGM^ /c and quadrupole moment qG^M^ jc^ 
{q < for oblate bodies) we can use post-Newtonian (PN) 
expansions (see, e.g., Weinberg 1972) to obtain approxi- 
mate expressions for the periastron advance and nodal pre- 
cession frequencies. Thus, for a nearly equatorial orbit of 
semilatus rectum I = 2r^T^I{r^ — rp), eccentricity e = 
(T'a — fp)/(T'a + 7'p), and the Newtonian and post-Newtonian 
orbital (Kepler) frequencies 

1 /l-^^^^''^ 



Vy. 



1 

2^ 

i^K.N 



I 



(1-36^)1 + 0^ ^ 



VP/2J 



(41) 



the rate of periastron advance is 



= I^K 



^1 + 2(1+6)^ 



^4j 



P/2 



3ql_ 
2 P 



+ 



[l^/^J 



(42) 



where O (-572) indicates terms of asymptotic order at least 

as high as l~"''^ . In equation ([42), the first two terms are the 
I''* and 2"'* PN corrections present even for a spherical star, 
the third term is due to the gravitomagnetic force, while the 
fourth term is the (Newtonian) effect of the star's rotation- 
induced deformation. 

In the lowest-order PN approximation, the precession of 
the plane of an orbit around the star's spin axis is interpreted 
as being caused by the gravitomagnetic force. Together with 
the purely Newtonian quadrupole precession, the nodal pre- 
cession rate is 

1 . 3q 1 „ / 1 



UK 



2j 



p/2 



-If-( 



15/2 



-J 



(43) 
^ [cf. equations (|37 



While for a Kerr black hole q 
and (By)], the ratio —q/j^ is considerably larger for neutron 
stars (Laarakkers & Poisson 1998); it ranges from ~ 2 to 
more than 10, depending on the mass and the equation of 
state. 



3 BOUND ORBITS AROUND NEUTRON 
STARS: COMPARISON WITH KERR 
BLACK HOLE ORBITS 

Equations of state of different stiffness give rise to a wide 
range of magnitude of the model neutron stars' rotation- 
induced deformation, and, consequently, significant differ- 
ences in the spacetime geometries around them. It is our 
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C 


uu 




L 






Re 

(km) 


(km) 


Re 

(km) 


(km) 


Re 
(km) 


^isco 

(km) 


3 


= 0.1 
= 0.2 


4.08 
(10.8) 

4.17 
(11.1) 


5.69 

(15.1) 

5.41 

(14.4) 


4.13 

(11.0) 

4.19 
(11.1) 


5.69 

(15.1) 

5.42 
(14.4) 


5.69 
(15.1) 

5.76 
(15.3) 


5.72 
(15.2) 






u^ 


9 


i/g 


g 


u^ 


9 


J 


= 0.1 
= 0.2 


271 
525 


-0.024 
-0.090 


249 
490 


-0.026 
-0.099 


148 
292 


-0.056 
-0.209 



Table 1. Stellar equatorial (circumference) radii Re, ISCO radii, 
spin frequencies u^ [Hz] and quadrupole moments q in units 
G'^M^/c'^ for NS models of mass M = 1.8 and indicated an- 
gular momenta j. The radii are given both in units GM/c^ and 
in kilometres. 



main purpose in this Section to describe and quantify how 
the differing geometries imprint themselves on the poten- 
tially measurable frequencies of bound orbits around the 
stars, and to compare these frequences with those of orbits 
around spinning (Kerr) black holes. 

For neutron star models of given mass M and angu- 
lar momentum jGM^ /c, the very stiff equations of state 
(e.g., EOS L as opposed to the 'softer' UU or yet 'softer' 
C) yield more extended neutron star models that undergo 
greater centrifugal deformation. This leads to higher mass 
quadrupole moments qG^M^ /c^ while requiring lower spin 
frequencies v^, as can be seen in Table 1. 

As Laarakkers & Poisson (1998, LP) have shown, 
the magnitude of q grows nearly quadratically, q ~ 
— o(EOS, M)j^, over a broad range of j. While the constant 
a is much smaller for EOS C (a ~ 2.4 for j <, 0.2; LP compute 
a ~ 2.7 for a much greater range of j) or EOS UU (a ~ 2.6) 
than for the stiff EOS L (a ~ 5.2; LP find a ~ 4.9), it still 
exceeds by a large amount the (exact) value a — 1 for the 
Kerr black hole. 

Based on the low-order PN expressions (tt2|) and (ksf), 
one would expect a general trend toward higher v^j^ (i.e., 
lower Vr) and lower fj^p at given r (all radii used in this sec- 
tion are the coordinate- independent circumferential radii) as 
—q is increased. This expectation is indeed borne out by our 
computation of circular orbital frequencies as can be seen, 
e.g., in Fig. 3 for a M = 1.8 family of NS models and black 
holes. Of course, a PN expansion cannot account for the 
strong-field relativistic effect of the vanishing of the radial 
epicyclic frequency at the radius risco of the innermost sta- 
ble circular orbit; the large magnitude of q for a stiff EOS 
simply causes a faster drop of z/^ as r is reduced past the 
point of maximum Pr- Quite generally, riaco moves inward 
for higher i, an effect familiar from the case of the Kerr 
spacetime (Ey). 

The vertical bars in Fig. 3 mark the ranges of the maxi- 
mal values of Ur (the top panel) and v^^-p (the middle panel) 
for the indicated values of j. For each j we plot the frequen- 
cies for the three EOS, C, UU and L (all solid lines), and 
the Kerr black hole (the dotted line) of the same Ad and j. 
In all cases the curves for EOS C and UU almost coincide. 




6 8 10 

r [GM/c'] 

Figure 3. The radial, nodal precession, and Kepler frequencies as 
functions of the circumference radius for circular orbits around a 
M = 1.8 neutron stars (solid lines; EOS C, UU,, and L) and Kerr 
black holes (dotted lines). The values of the angular momentum 
parameter j are indicated by the labels of vertical bars. 

Whereas the drop in Vr sd, r < 10 due to the star's spin- 
induced deformation is relatively modest [\ Aurl / i^r ^^5%; 
notice, however the considerable numerical error in u,. for the 
stiff EOS L, due to the necessity of computing the second 
derivative (nsl) from the tabulated metric], the reduction 
in Vf^pir) relative to the Kerr values is more pronounced, 
amounting at j — 0.2 to up to ~ 20% for EOS UU or even 
~ 40% for EOS L. For lower j {j < 0.1), the difference be- 
tween the frequences for different EOS is less pronounced, 
and it virtually disappears at j ~ 0.05. 

Of the three EOS, the stiff EOS L makes the star at 
j = 0.2 so extended at the equator. Re = 5.8GM/c^ = 
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D. Markovic 
0.2 

0.1 



-0.1 



-0.2 



U q= 



-0.21 



8 10 

r [GM/c''] 



20 



Figure 4. Relative deviation, (i^p^ — !^pa k )/'^pa k > with respect 
to the Kerr values (for the same M and j) of the periastron ad- 
vance frequencies (solid line) and the post-Newtonian values (dot- 
ted line; see equation l43) for circular orbits around an M = 1.8, 
j = 0.2, EOS L (g = -0.21) neutron star. 




-0.1 



-0.2 



-0.3 



-0.4 



6 8 10 20 

r [GM/c''] 

Figure 5. Relative deviations, (i^j^p — ^^np k)/'^np k 
spect to the Kerr values (for the same M and j) of the nodal 
precession frequencies (solid lines) and the post-Newtonian val- 
ues (dotted lines) for circular orbits around M = 1.8, j = 0.2 
neutron stars. The EOS used are UU (yielding q = —0.10) for the 
upper solid line, and L {q = —0.21) for the lower solid line. The 
post-Newtonian values are obtained from equation (Usb using the 
values of M, j and q of the corresponding neutron star models. 



with re- 




4 6 8 10 20 

r, [GM/c^] 

Figure 6. The radial, nodal precession and Kepler frequencies as 
functions of the apastron radius for eccentric orbits of common 
periastron rp = 5.88 around M = 1.8 neutron stars (solid lines; 
EOS C, UU and L) and Kerr black holes (dotted lines). The values 
of the angular momentum parameter j are indicated by the labels 
of the vertical bars. 



15.4 km, that there is no innermost stable circular orbit (see 
Table 1); the corresponding curve Vrir) thus never reaches 
the terminal point i/r ~ 0, in contrast with all the other 
cases shown in Fig. 3. 

The simple low-order PN expressions (Ud) for i/p^ and 
( [431) for Uj^p approximate rather accurately the actual fre- 
quencies of circular orbits of radii r ^ 10, as can be seen from 
Figs. 4 and 5. As pointed out above, close to riaco (i-e., for 
r < 8), the growth of the post-Newtonian Vpj^ with decreas- 
ing r cannot model accurately the much faster growth of the 
relativistic //p^ as it approaches uk- In this range, the Kerr 



values for the same M and j furnish a considerably better 
approximation, as is evident in Fig. 4. 

For the nodal precession, on the other hand, the large 
quadrupole term makes at all r the low-order expansion mM 
a much better approximation than the Kerr i^j^p{r) (see 
Fig. 5). Notice, however, that for stiff EOS and for r<,8 
(the lower curves of Fig. 5) the expansion (US) significantly 
overestimates the actual v^p which is, e.g., at the surface of 
the star, less than 60% of the corresponding Kerr value. 

Turning now to eccentric orbits, one can distinguish be- 
tween two broad classes, characterised by periastra Vp that 
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0.4 
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0.008 
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Figure 7. Ratios of orbital frequencies for M = 1.8, j = 0.2 and 
rp = 5.88 for EOS L and UU (solid lines) and for the Kerr black 
hole (dotted lines) of the same M and j. PN values for q = —0.21 
(EOS L) and q = -0.10 (EOS UU) are given by dashed lines. 
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5 10 50 50 100 

r\ [GM/c^] r; [GM/c^] 

Figure 8. Relative deviations, (i^p^ — '^pa k)/'^pa k> with re- 
spect to the Kerr values (for the same M and j) of the peri- 
astron advance frequencies (solid lines) and the post-Newtonian 
values (dotted lines; see equation [42| ) for eccentric orbits around 
an M = 1.8, j = 0.2, EOS L (g = -0.21) neutron star. Left 
panel: eccentric orbits of common periastron rp = 5.88. Right 
panel: eccentric orbits of common periastron rp = 15.0. 

are either larger or smaller than riaco- Within the first class, 
^p > ^isco (so that ra.min = ''p), the plots of the orbital fre- 
quences as functions of r^ for a family of orbits of common 
periastron rp resemble (an example is shown in Fig. 6) the 
frequency vs. radius plots for circular orbits shown above in 
Fig. 3. The obvious difference is that the magnitudes of the 
negative slopes diz/dra. are smaller than those of the slopes 
du/dr for circular orbits: even if ra is large, the eccentric or- 




5 10 50 50 100 

r, [GM/c^] r^ [GM/c^] 

Figure 9. Relative deviations, {v-^-p — v^^ k)/'^np ki ™ith re- 
spect to the Kerr values (for the same M and j) of the nodal 
precession frequencies (solid lines) and the post-Newtonian val- 
ues (dotted lines; see equation |43| ) for eccentric orbits around 
M = 1.8, j = 0.2 neutron stars. The EOS used are UU (yielding 
q = —0.10) for the upper solid and dotted lines, and L (q = —0.21) 
for the lower lines. Left panel: eccentric orbits of common perias- 
tron rp = 5.88. Right panel: eccentric orbits of common periastron 
rp = 30.0. 



bits have access to small radii where relativistic effects, like 
frame dragging, take plac e. This can be seen from the PN 
expansions (|42| ) and (|43|), in which the effect of increased 
periastron radius is 'buffered' by ra's appearing in the ex- 
pressions only through I, where / ~ 2rp for ra ^ rp [the 
dependence on e in (tt2|) is weak]. For large ra we therefore 
expect the ratios v^f^/uK and v^^/vk to asymptote to con- 
stant values, as indeed is the case (see Fig. 7). 

Although the PN expressions (W2|) and ( |43| ) reproduce 
qualitatively the profiles i/p^(ra) and v^-pij'^) obtained nu- 
merically from NS models, they tend to give a somewhat 
steeper fall-off at large ra, thus falling short of the 'exact' 
values by up to ~ 30%. This is evident in Fig. 7 and is dis- 
played in greater detail in Figs. 8 and 9. It is a consequence 
of the extreme slowness of relativistic radial motion in the 
strong-filed vicinity of risco: since relativistic effects (e.g., 
frame dragging) weaken rapidly at larger radii, this gives 
them more time to build up and the integrated effects thus 
exceed the PN estimates. 

As in the case of circular orbits, the periastron advance 
frequency J^pA('"a) is modeled much better for all ra by the 
Kerr black hole values (for the same M and j; see Fig. 8) 
than by the PN expansion ( |42[ ) . As can be seen in the right 
panel of Fig. 8, only for rp > 10 does the PN i-'pA(ra) 
approach closely the values obtained from the numerically 
computed NS spacetimes. 

Regarding the nodal precession, for moderately eccen- 
tric orbits, ra^2rp, the PN expansion (Esl) is superior to 
the Kerr values at modeling v-^-pirs,) obtained from the nu- 
merical NS spacetimes (see Fig. 9). The better accuracy of 
the PN frequences is more prominent for stiff EOS, in a 
straightforward extension of the trend we observed for cir- 
cular orbits. Again, for larger eccentricities, the numerical 
spacetimes' t'Np(^a) inevitably split off from the PN curves 
and shift toward the Kerr values. As we increase rp, the 
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15 20 25 

[GM/c^] 



Figure 10. The radial, nodal precession and Kepler frequencies 
as functions of apastron radii of highly eccentric orbits of common 
periastron Tp = 4.30 around M = 1.8 neutron stars (solid lines; 
EOS C and UU) and Kerr black holes (dotted lines). The values 
of the angular momentum parameter j are indicated in the top 
panel. 



range of eccentricity e over which the PN expansion is rela- 
tively accurate increases; compare the two panels of Fig. 9. 
For the second class of eccentric orbits with Tp < risco, 
we observe (see Fig. 10 for EOS C and UU) that at constant 
Tp the minimum apastron radius, ra,min(7'p), shrinks rapidly 
with increasing j, following the related reduction of risco (for 
both EOS, risco = 5.84, 5.69 and 5.42 for j = 0.05, 0.1 and 
0.2, respectively). Compared with the Kerr black hole, the 
main effect of the larger quadrupole moment of the NS is a 
larger ra,min(rp) at given M and j as well as lower t'r('"a) and 
higher i^p^(ra). Notice, however, that as ra approaches ra,min 
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Figure 11. EOS UU; frequencies and apastra as functions of 
J^K for families of eccentric orbits around a NS of M = 1.7 and 
i/^ = 100 Hz (Re = 4.39 = 11.0km, Hsco = 5.87 = 14.7km). The 
curves are labeled by the values of periastron in units GMjc? . 
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Figure 12. EOS L; frequencies and apastra for families of eccen- 
tric orbits around a NS of M = 2 and v^ = 363 Hz (i?c = 5.23 = 
15.4km, risco = 5.44 = 16.1km). 
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from above, the functions !/K(ra) and i^^p{ra) bend upward 
abruptly, while z^r(''a) drops until it reaches the 'terminal' 

value i^r(r-a,min)- 

The sharp change of i/k near ra.min leads to a 'flat' de- 
pendence Vr{i'K) for a family of orbits of a given r^ < risco; 
see Fig. 11 for a relatively slow spinning, v^ — 100 Hz, EOS 
UU neutron star. As rp is increased, the functions Vr(i^K) 
at first grow steeper, with increasingly negative slopes at 
large vk (small ra), until at Vp — risco, i^r drops to zero for 
fK = i-'isco- A further growth of rp again brings about flatter 
functions i'r{i'K) but at considerably larger z^,.- 

In Fig. 12 we show the functions Ur{uK) for eccentric 
orbits around a. M = 2, v^ = 363 Hz, EOS L neutron star. 
The stiff equation of state makes the star large (-Re = 5.3 = 
15.6 km) so that very low rp orbits do not exist and the flat, 
low Vr curves Uriya) seen in Fig. 11 are absent from Fig. 12. 
(This is why the plots of frequencies for the low rp = 4.3 
families of orbits in Fig. 10 could be shown only for the 
moderately stiff EOS C and UU.) As a consequence, the set 
of frequencies of eccentric orbits (for given M and i^^) is a 
much more limited subset of the 2-dimensional v^-Vr plane 
for the stiff EOS L than for the softer equations of state. 



4 CONCLUSIONS 

If shown to indeed be a manifestation of the orbital mo- 
tion around neutron stars, the quasi-periodic oscillations of 
the low mass X-ray binaries would furnish an unprecedented 
probe of the strongly curved spacetime near the compact ob- 
jects and allow inference of their masses M and spin rates 
i/g. In particular, highly eccentric orbits, which near their 
periastra penetrate deeper than the innermost stable circu- 
lar orbit, would put tight constraints on the sizes of neutron 
stars, and thus on the equation of state. 

The high accuracy (e.g., better than ~ 1% for Sco X-1), 
with which the frequencies of the kHz QPO frequencies have 
been measured, makes it important to use the values ob- 
tained by integrating test-particle motion in fully relativis- 
tic numerical spacetimes that can be computed by spinning 
neutron star codes. For j > 0.1, and especially for orbits of 
considerable eccentricity, the post-Newtonian (FN) expan- 
sions or approximation with Kerr metric orbital frequencies 
are not sufBciently accurate. 

In the case of the periastron advance frequency u-pj^ , the 
Kerr results provide a better approximation to actual motion 
in the neutron star spacetimes than do the low-order FN ex- 
pressions. Although the Kerr z/p^ j^ is always lower than the 
value fpj^ for the same periastron and apastron (circumfer- 
ential) radii in a neutron-star spacetime, it never falls below 
Vpj^ by more than 10%, even for the highest angular mo- 
menta investigated here (j « 0.2) and for the stiffest EOS 
(L). The general-relativistic decrease of Vr in the vicinity of 
the innermost stable orbit (ISCO) makes the FN approxi- 
mations for Vpj^ inaccurate even for moderately (r^/rp < 2) 
eccentric orbits if rp < 7 (for highly eccentric orbits, FN is 
typically a poor approximation). 

The nodal precession frequency !/„p of orbits inclined 
relative to the equatorial plane of a neutron star is domi- 
nated by the prograde gravitomagnetic 'torque,' with a large 
retrograde (for prograde orbits) contribution from the clas- 
sical quadrupole torque. In the cases considered in this pa- 



per, the latter torque can amount at small orbital radii from 
about 10% to up to about 50% of the magnitude of the grav- 
itomagnetic torque. The spin-induced deformation (oblate- 
ness) of the neutron stars renders, therefore, the low-FN z/„p 
far more adequate (for orbits of moderate eccentricity) than 
the values obtained from Kerr black holes with the same M 
and j. However, for orbits of very small size (ra ^8), even 
the low-order FN expressions overestimate v^^-p by up to 30% 
in the case of j ~ 0.2, stiff EOS L models. 

Very stiff equations of state (e.g., L) have an impor- 
tant effect on the high-frequency portion of the set of all 
available frequencies z/r, I'r of eccentric orbits: due to the 
large stellar equatorial radii, orbits with very low perias- 
tra do not exist. The consequent absence of low values of 
Ur at t'K,isco/2 < !^ < z^K.isco, for given M and u^, restricts 
severely the available subset of the vy^-Vr plane. For this 
reason, very stiff equations of state do not permit successful 
modeling of the neutron-star LMXBs' kHz QPO frequencies 
in terms of the frequencies of eccentric orbits (see Markovic 
&Lamb 2000). 
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